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A fluid flow in a simple dense liquid, passing an obstacle in a two-dimensional thin film geometry, 
is simulated by Molecular Dynamics (MD) computer simulation and compared to results of Lattice 
Boltzmann (LB) simulations. By the appropriate mapping of length and time units from LB to 
MD, the velocity field as obtained from MD is quantitatively reproduced by LB. The implications 
of this finding for prospective LB-MD multiscale applications are discussed. 
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Micro and nano-^hydrodynamic flows arc playing an 
increasing role for many applications in material sci- 
ence, chemistry, and biology P|. To date, the leading 
tool for the numerical investigation of nano and micro - 
hydrodynamic flows is provided by molecular dynamics 
(MD) simulations P|. In principle, MD yields a correct 
description of fluids on microscopic and hydrodynamic 
scales. But typical length and time scales that can be 
covered by MD simulations of dense liquids are of the 
order of a few tens of nanometers and about a few hun- 
dred nanoseconds, respectively. As a result, the quest for 
cheaper — and yet physically realistic alternatives — is 
a relentless one. Mesoscopic models, and notably those 
arising from kinetic theory, are natural candidates to fill 
this gap because they operate precisely at intermediate 
scales between the atomistic and continuum levels. Very 
recently, the so-called Lattice-Boltzmann (LB) method, 
where a kinetic equation is solved on a lattice, has proven 
successful to model characteristic features of microscopic 
flows, such as the occurrence of slip boundary conditions 
near solid walls Q . 

However, whether LB is also able to reproduce complex 
nanoscopic flows on a quantitative level, in particular for 
a dense liquid, remains an open issue. Although, pioneer- 
ing work of various authors has shown that hydrodynam- 
ics holds nearly down to the molecular scale in simple 
dense liquids (see, e.g. 01), it is not obvious whether in 
the vicinity of wall-fluid or obstacle-fluid boundaries at 
least pair interactions on an atomistic level have to be 
taken into account to yield a quantitative description of 
fluid flows. Although LB is a particle-based method, 
it cannot resolve the short ranged structural order of a 
dense liquid since it describes a structureless lattice gas 
with no many-body interactions. Thus, it is not clear 
whether one can match the LB method with MD simula- 
tions of complex flow structures in dense liquids (e.g. co- 
herent nanostructures). 

In this paper, we aim to establish a quantitative map- 
ping between the resolution requirements of a LB ver- 
sus MD simulation of a non-trivial nano-hydrodynamic 



flow. Far from being a purely academic exercise, this is 
a primary issue to set a solid stage for future multiscale 
applications combining LB with atomistic methods. To 
this end, we have performed MD simulations of a two- 
dimensional Lennard- Jones fluid confined between walls 
that passes a thin plate-like obstacle while subject to a 
gravitational force field. This system is then modeled 
by the LB method. We demonstrate that, after the ap- 
propriate mapping of length and time units, there is a 
quantitative agreement between the flow fields computed 
with MD and LB. 

First, we give a brief introduction to the LB method 
and the LB model used in this work. The LB method 
as a mesoscopic simulation tool has met with significant 
success in the last decade for the simulation of many 
complex flows 0. We consider the LB equation in its 
original matrix form [^: 

/,(x + Atc„t + At)-/,(x,t) = 

-^%At [/,(x,i)-/;(x,t)] (1) 

3 

where Ai is the time unit and fi{'x.,t) = /(x, v = Ci,t), 
i = 0, . . . n, is the probability of flnding a particle at 
lattice site x at time t, moving along the lattice direc- 
tion defined by the discrete speed c^. The left-hand 
side of Eq. represents the molecular free-streaming, 
whereas the right-hand side represents molecular colli- 
sions via a multiple-time relaxation towards local equi- 
librium /° (a local Maxwellian expanded to second order 
in the fiuid speed). The leading non-zero eigenvalue of 
the collision matrix fixes the fluid kinematic viscos- 
ity as ly ~ c^{l/uj — 1/2) (in lattice units At = Ax = 1), 
where Cg is the sound-speed of the lattice fluid, l/VS in 
the present work. In order to recover fluid dynamic be- 
haviour at macroscopic scales, the set of discrete speeds 
must be chosen in such a way as to conserve mass and 
momentum at each lattice site. Once these conserva- 
tion laws are fulfilled, the fiuid density p = J2i hi ^'^d 
speed u = J2ifi^i/P can be shown to evolve according 
to the quasi-incompressible Navier-Stokes equations of 
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fluid-dynamics . 

In the bulk flow, LB is essentially an efficient Navier- 
Stokes solver in disguise. The major advantages over a 
purely hydrodynamic description are: i) the fluid pres- 
sure is locally available site-by-site, with no need of solv- 
ing a computationally demanding Poisson problem, ii) 
momentum diffusivity is not represented by second order 
spatial derivatives, but it emerges instead from the first- 
order LB relaxation-propagation dynamics. As a result, 
the time-step scales only linearly (rather than quadrati- 
cally) with the mesh size, which is an important plus for 
down-coupling to atomistic methods, iii) highly irregu- 
lar boundaries can be handled with ease because particles 
move along straight trajectories. This contrasts with the 
hydrodynamic representation, in which the fluid momen- 
tum is transported along complex space-time dependent 
trajectories defined by the flow velocity. 

At the fluid-solid interface, we adopt the no- slip 
boundary condition, u = 0. This is imposed by the stan- 
dard bounce-back procedure, i.e. by reflecting the out- 
going populations back into the fluid domain along their 
specular direction. With reference to particles propagat- 
ing south-east(\) from a north-wall boundary placed 
a.t z = H + 1/2, the bounce-back rule simply reads as 
f\{x,y,H+l) = f^{x + l,y,H) where the lattice spac- 
ing is made unity for convenience. This corresponds to 
a stylized two-body hard-sphere repulsion, with interac- 
tion range equal to y/2 lattice units. 

In this work, we use a standard nine-speed D2Q9 
model [3 to simulate a two-dimensional channel flow of 
size L and H along the streamwise and cross-flow direc- 
tions, respectively. At a distance L/A from the inlet, a 
thin flat plate of height h is placed, perpendicular to the 
direction of the main flow. The fluid is driven by a con- 
stant body force fx in x-direction which, in the absence 
of the vertical plate, would produce a parabolic Poiseuille 
profile of central speed U = fxH^ / {^v), v being the kine- 
matic viscosity of the fluid. 

A similar two-dimensional channel flow of a fluid is 
considered in the MD simulations. As a model for the 
interactions between the particles, a Lennard-Jones po- 
tential is used that is truncated at its minimum and then 
shifted to zero. It has the following form: 



V(r) 




7)" -(7)^ 



for r < 2I/V 
otherwise, 



(2) 



where r denotes the distance between two particles. The 
units are chosen such that e, cr, the Boltzmann constant 
/cb and the mass of the particles m are set to one. The 
particles are confined into a rectangular box of size Lx H 
with L = 200tT and H = 106a. Moreover, an obstacle is 
placed along a line at a; = L/4. This obstacle consists 
of 41 fixed particles between y = 33.0a and y = 73.0(7 
with an equidistant spacing of Icr. Thus, the obstacle 
has an effective width of about la and a height h = 



AOa. For the density we chose n = N/{LH) = 0.8cr~^ 
corresponding to = 16879 particles in the simulation 
box. The equations of motions were integrated using 
the velocity form of the Verlet algorithm with a time 
step of dtmd = O.OlTmd with = \/ ma^ / (48e). In 
order to keep the temperature constant, a Nose-Hoover 
thermostat was applied |^. 

The MD simulations were all done at temperature 
T ~ 5.3e/fcB. Ten independent runs were performed in 
order to obtain reasonable statistics. First, the systems 
were equilibrated for 30000 time steps. Then, walls were 
introduced in the system by giving all the particles at 
X < 3.0(7 and at a; > 103.0(7 zero velocity. The choice 
of these rough walls provides the absence of any layer- 
ing effects of the fluid near the walls. The system with 
walls, i.e. with the immobilized particles, was then equi- 
librated for another 20000 time steps, followed by runs 
over two million time steps with a gravitational force fleld 
perpendicular to the walls {x direction) with magnitude 
fx = 2x 10~^ at each particle. Note that in the MD with 
external force fleld the Nose-Hoover thermostat was only 
applied in y-direction, i.e. perpendicular to the flow field. 
Results were collected after 2 x 10^ time steps when the 
steady state was clearly reached. Then, every 1000 steps 
a configuration with positions and velocities of the par- 
ticles was stored from which all the quantities that are 
presented in the following were computed. In total the 
average was over 18000 configurations. 

In addition to the latter runs, also MD simulations 
were performed for a system without obstacle. Here, the 
aim was to determine the kinematic viscosity from the 
Poiseuille proflle that forms at steady state. As a result, 
we obtained the value — 0.5 for the kinematic viscosity 
(in MD units). Thus, with U{ = 0.07 the average flow 
velocity in the simulations with obstacle, the Reynolds 
number is Re = Ufh/v « 5.6, above the critical value for 
the onset of coherent structures. 

In order to compare the results of the MD with those 
of LB, a conversion of MD units into LB units is required. 
Space and velocity (time) conversion proceed as follows: 
L = LibAx = Lnidcr, U = UihAx/At = L^^a/T^A- This 
yields Lib = imdcr/Aa; and t/ib = t^md^-^- The con- 
version of the kinematic viscosity, central to the definition 
of the Reynolds number, is given by v\h = t^md(3j)^:^- 

Reducing Aa; to values of the order of a means that the 
LB simulation would resolve the structure of the fluid, if 
this happened to be included into the lattice kinetic equa- 
tion. This is not the case here. The interesting question, 
though, is whether the absence of these structural effects 
does hamper the correct reproduction of bulk-flow fea- 
tures. This is very similar, although on totally different 
scales, to the issue of subgrid scale modeling of turbulent 
flows: do the unresolved scales spoil the physics of the 
resolved ones? 

In order to map LB units onto MD units, wc proceeded 
as follows: First, wc identified the LB mesh-spacing with 
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FIG. 1: Steady-state streamwise velocity profiles u{y) —< 
Vx{y) > as a function of the crossflow coordinate y at x = 
L/4 + 5(7 (i.e., 5a behind the obstacle) and x = L/2 as ob- 
tained from the MD as indicated. The solid lines are the cor- 
responding results from the LB with resolution Ax = 0.25cr. 
Both the MD and LB curves are averaged along the x direc- 
tion over a region of ±3(t. The inset shows a magnification of 
the region around the minimum for the profile at a; = L/4-f 5(t 
(here, the dashed line is the LB result without the averaging 
along the x direction. 



a fraction or a multiple of a. Then, the time conversion 
factor was determined from the kinematic viscosity, such 
that the same Reynolds number. Re ~ 5.6, was yielded 
in LB and MD. In this way, LB simulations were done 
for Ax = 0.25 tr, 0.5 cr, LOtr, and 2. Oct which respectively 
corresponds to the values i>\i, = L92, 0.96, 0.48, and 0.24. 
The LB simulation ran over 10000 steps, which was found 
sufficient to keep time changes of the overall velocity pro- 
file at least within third digit accuracy. The CPU time re- 
quired for the LB simulations varied between about 1 min 
to 45 min on a Pentium 4 with 2.8 GB clock rate, depend- 
ing on the chosen resolution. This has to be compared 
to the total computational load for the MD simulation, 
which was about 1 week on 10 Pentium 4 processors. 

In Fig.^ we show the steady-state streamwise velocity 
profile u{y) =< Vx{y) > as a function of the crossflow 
coordinate y at a: = L/4 -I- 5ct and x = for the 

LB with Ax — 0.25ct, as compared to the corresponding 
profile obtained with the MD simulations. Here, both 
MD and LB data are averaged along the x direction over 
a region of ±3ct. This leads to a smoothening of the MD 
data, but, as the LB results show, has only minor effects 
on the profiles (this can be infered from the dashed line 
in the inset of Fig. ^ which shows the LB data without 
averaging along the x direction). From the data, very 
good agreement is observed between MD and LB results. 
Of particular interest is the region around the minimum 
in the curve for x = L/4 -I- 5ct, where the velocity field 
< Vx > becomes negative. These negative velocities are 
due to the formation of vortices behind the obstacle. As 
the inset of Fig. ^ shows, even this region of negative 
velocities is well reproduced by the LB simulation. 



FIG. 2: Steady-state streamwise velocity profile u(y) —< 
Vx(y) > as a function of the crossflow coordinate y a.t x = 
L/A + 5a and x = L/2 for the LB simulations with different 
choices of Ax as indicated. As in Fig. the inset shows a 
magniflcation of the region around the minimum in the curves 
for X = L/4 + 5a. 
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FIG. 3: Steady-state streamwise velocity profile u{y) —< 
Vx(y) > as a function of x at y = H/2 and y = H/2 ± 25a. 
The symbols show the MD data as indicated. The solid lines 
are the corresponding LB data with Aa; = 0.25cr. 

In Fig. [5] the same quantities as in Fig. ^ are shown 
for the LB runs with different choices of Ax. From this 
figure, it is apparent that lack of resolution of the atomic 
scale CT generates significant departures of the LB results 
from the MD data. This is an informative result, for it 
implies that microscopic length scales, although absent 
from the LB description, must nonetheless be resolved if 
the flow proflle away from the walls is to be quantita- 
tively captured by the coarse-grained simulation. This 
is especially evident in the region where vortices form, 
i.e. just behind the obstacle. In this region, a quantita- 
tive agreement with the MD requires a resolution as high 
as Ax = 0.25ct (see the inset of Fig.|2J. 

In Fig. 121 the streaming velocity < w^, > in the direc- 
tion parallel to the flow is considered for y = H/2 and 
y = H/2 -\- 25ct. In this case LB and MD data were av- 
eraged over a region of ±2ct. Note that in the MD case, 
also the data &i y ~ H/2 ~ 25ct were used for the average 
of the profile at y = H/2 + 25ct, since the profile is ex- 
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FIG. 4: Upper panel: The color map shows the magnitude of 
the velocity field, |u(2;, = < Vx >^ + < Vy >^ at steady 
state as obtained from the MD simulation. Lower panel: The 
same result as obtained from the LB simulation with Aa; — 
0.25a is shown. 

pected to be symmetrical with respect to the central line 
at y = H/2. The comparison between LB and MD shows 
again a very good agreement. This is not a foregone re- 
sult, since there are non-trivial features in the vicinity 
of the obstacle, i.e. around x = 50a, such as the shoul- 
der at a; = 55(t in the curve for y = H/2 + 25a. Fig. 0| 
displays the magnitude of the velocity field, |'u(a;,y)| = 
y/ < Vx >^ + < Vy as obtained from MD (upper plot) 
and LB (lower plot). Here, it is illustrated that the whole 
velocity field is quantitatively reproduced. 

In conclusion, we have demonstrated an example of a 
complex nanoscopic fluid flow, in which a spatial hand- 
shaking between LB and MD is possible. The mapping 
of LB onto MD requires a conversion of time and length 
units from LB to MD and one has to choose an appropri- 
ate grid resolution for the LB fluid. For the dense fluid 
considered in the MD simulation of this work, a grid res- 
olution Aa; = 0.25(7 is required to reproduce the flow 
features quantitatively. Remarkably, non-hydrodynamic 
(finite-Knudsen) effects appear to be silent, at least at 
steady state. This is likely to be a benefit of the matrix 
formulation of the collision operator 0, j as compared 
to the more popular single-time relaxation form. 

The present results indicate that there appears to be 
a sound ground for prospective multiscale applications 
based on the combined use of (multigrid) LB [13 with 
MD. For instance, by using multigrid LB with, say, 6 
levels of resolution, one could couple LB with MD at the 
finest scale Aa;finc ^ a, typically near the boundaries, 
and then progressively increase the LB mesh-size so as 



to reach a hundred-fold larger mesh spacing Axbuik = 
2^Ax{ine in the bulk flow. On-the-fly time-coupling is 
obviously more demanding. In fact, our simulations in- 
dicate that even with the finest resolution Ax = a/4, LB 
is about a factor 2000 faster than MD (see above). The 
most direct strategy for 'on-the-fly' LB-MD coupling is 
to apply LB everywhere in the fluid domain and leave 
MD only in small portions, typically in a ratio f : f 000 
to the global domain. Moreover, this factor fOOO gap 
could be partially bridged by resorting to very recent 
time-adaptive LB procedures This indicates that, 
although very demanding, even time-coupling between 
LB and MD may become soon feasible for the numerical 
investigation of complex nanoflows. 

Of course, several open issues remain. Our compar- 
isons were done at steady state and it is not clear to 
what extent transient states are correctly described by 
LB. Moreover, it will be interesting to see whether the 
LB method still yields a quantitative description in the 
case of slip-flow and other complex phenomena at the 
fluid-solid interface. These issues make interesting sub- 
jects for forthcoming studies. 
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